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bubble populations in a moving fluid. The bubble distribution, as a function 
of position, velocity, radius, and time, satisfies a Boltzmann type transport 
equation that is derived and then formally solved by the method of character- 
istics. Before this new analytical tool can be applied, properties of the 
medium must be specified and a bubble dynamics model must be chosen. General 
expressions are written for bubble acceleration and radius change rate, and 
known models of bubble gas diffusion and drag are summarized for gas bubbles 
in liquids. Application to the upper ocean is discussed and then illustrated 
with seme sample calculations. Explicit solutions are written for a steady- 
state, one dimensional ocean with distributed sources; the results clarify 
relations between observed bubble populations, proposed bubble source mechanisms, 
and knewn models of single bubble dynamics. 
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I. 



INTRODUCTION 



Ihe solutions of many experimental and engineering problems depend a 
knowledge of bubble distributions in a fluid medixim. For instance, gas bubbles 
near the surface of the ocean are important in a variety of subjects that 
include underwater sound propagation (Shulkin 1968 & 1969) , meteorology 
(Blanchard & Woodcock. 1957) , sea surface chemistry (Sutcliffe, et al 1963), 
cavitation (Rdx and Herzfeld 1954) , and air-sea gas exchange (Kanwisher 

1963) . Several efforts have been made to measure bubble densities in 
the ocean (Blanchard & Woodcock 1957; Glotov, et al 1962; McCartney and Bary 
1965) , the most recent made by Medwin (1970) . To infer near surface bubble 
populations from in-situ measurements of acoustic attenuation, Medwin exploited 
the fact that a gas bubble's scattering and absorption cross sections for sound 
at the bubble resonant frequency can typically be 1000 times its gecmetrical 
cross section. For this reason large bubble densities near the sea surface can 
significantly affect xjnderwater sound transmission. 

Many corplex factors affect bubble distributions in a fluid such as the 
ocean, and some of these phencmena have been carefully investigated. For 
example, LeBlcnd (1969) examined gas diffusion from an ascending bubble, and 
Levich (1962, Ch. 8) is an excellent reference on single bubble dynamics. The 
purpose of this paper is to develop and to demonstrate new analytical tools 
for investigating relationships between observed bubble distributions and the 
phenomena vdiich create, transport, and eliminate bubbles. 

In sections two and three the formalism of transport theory is adapted to 
the general problem of describing bubble populations in a moving fluid. The 
bubble distribution, as a function of position, time, velocii^, and radius, 
satisfies a Boltzmann-type transport equation that is derived. Using a 
Lagrangian viewpoint, the formal solution is expressed as a line integral of 
the volume source function along characteristic curves. In section four general 
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expressions are obtained for bx±ble acceleration and radius change rate and 
sane known nodels of gas diffusion and bubble drag are introduced. Finally, 
in the context of gas bubble transport in the upper ocean, a simple model is 
chosen to illustrate the utility of transport theory. Calculated results 
clarify relationships between observed bubble distributions, proposed bubble 
source mechanisms, and known models of single bubble dynamics. 
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II. THE BUBBLE TRANSPORT EQUATION 



Consider the problem of describing the bubble distribution — as a function 
of position, r, velocity, v, size, £, and time,,t— in a moving fluid that may 
contain bubble soxirces and sinks. In general, a "bubble" might be considered 
any sinple, closed region containing a ’fluid somehow different from that of 
the transporting medium. Since exaiiples chosen in this paper refer primarily 
to spherical gas bubbles, i will represent the bubble radius. However, if the 
bubbles are not spherical, then £ might represent the radius of a spherical 
volxmie egucil to the actual bubble vol\jme, or one might wish to generalize the 
description to include more than one parameter for size. 

Suppose that the time evolution of the position, the velocity, 
and the radius of each bubble depends only upon the properties 
of the medium and the bubble characteristics and is statistically independent 
of the state of any other bubble. Then define the bubble distribution function, 
i|;(r,v, £,t)d3rd^vd£, as the mean number of bubbles at time t in the volume d^r 
about r, with velocity in dV about v and the radius in d£ about £.* Further- 
more, if bubbles interact with one another, suppose that no more than t3VO 
interact at any one time and that the collision time is short corpared to the 
time for ip to change appreciably. Then ip satisfies a Boltzmann-type transport 
equation whose solution describes the ensemble average behavior of the bubble 
population in terms of single bubble behavior. 

At this point it is convenient to define the seven-dimensional volume ele- 
ment d^T H d^rd^vd£ and to adopt a Lagrangian viewpoint. In the time interval 
(t,t+dt) the bubbles that were in d^x about (r,v,£) have moved to a new volume 



*E.g., in the Cartesian representation, position components lie in the intervals 

(x,x+dx) , (y,y4dy) , and (z,z4dz); velocity components lie in the intervals 

(VxfVx+dvx) , (Vy,v„4d.Vy) , and (v 2 ,V 2 +dv 2 ) ; and the radius is in the interval 
(£,£+d£). -1^ y -S' 
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d^T about (r',v',£'), where the primed quantities refer to the time t+dt and 
the unprimed quantities refer to the time t. This transport occurs by a change 
in position due to a velocity, v = dr/dt, by a change in velocity due to an 
acceleration, a = dv/dt, and by a change in radius due to a rate of change, 

V = di/dti 



r' = r + vdt + O(dt^) 

v' = V + adt + O(dt^) (1) 

V = I + vdt + 0(dt2) 

(Taylor's series about time t.) For gas bubbles in the ocean and elsevi^ere, 
the acceleration results frcsn several influences that include gravity, buoyancy, 
and entrainment of bubbles by the fluid, vdiile the rate of change of bubble 
radii is caused by hydrodynamic canpression, surface tension and gas diffusion. 
Expressions for a and v are obtained in section four vdiere single bi±>ble 
dynamics is discussed. 

Conservation of bubbles denands that the number of bubbles in the volume 
element at time t+dt, i|;{r' ,v' ,£' ,t+dt)d^x' , equals the number that start in the 
volume element at time t, tj;(r,v, Z,t)d^T, plus (minus) any bubbles that are 
introduced (lost) by distributed sources (sinks), S(r,v,£,t)d^xdt: 

i|)(r' ,v' ,2.' ,t+dt)d^x' - i|)(r,v, 2,t)d^x = S(r,v,£,t)d^xdt + O(dt^), (2) 

where S(r,v,2,t)d'^x is the rate at vdiich bubbles are introduced into d^x about 
(r,v, 2) distributed soxirces. To develop equation (2) it is sufficient to 
note that terms of O(dt^) will disappear in the limit dt -»■ o. Therefore we 
expand ij)' e i|)(r' ,v' ,£' ,t4dt) in a Taylor series about i|» = i|;(r,v,2,t) and keep 



8 



terms through 0 (dt) : 



ip' = Ijj + (,8i|)/8t)dt + (v-v^i|j)dt + (a.^^i/))dt + vOip/3Jl)dt + 0(dt2) (3) 

where V^(-) = (9(-)/9x,9(*)/ay,3(-)/3z) and V^(.) = (9 (-)/9v^, 3 (•)/9v^, 3 (OAv^) 

are Cartesian representations of the "del" operators in position space and 
veloci-tY space, respectively. 

To proceed, d^x ' must be expressed in terms of d^x : 

d^x' = Y(t4dt,t)d”^x (4a) 

where 

3 (r * ,v' ,Z') 

Y(t4dt,t) s H det ((q^qj^)) (4b) 

i,j — 1,2,«.«,7 



is the Jacobian with 

^1 = ^2 - ^4 - ^ ^ z ' ^7 - 

q^' = x', q2* = y', etc. 

The expression derived for y in the appendix can be expanded in a Taylor series 
about t, or equation. (1) can be used to ccitpute the determinant in (4b) directly. 
In either case the result is 

Y(t-hlt,t) = 1 + (V^-a)dt + (3v/3Jl)dt +0(dt^) (6) 

Substituting (3), (4) and (6) into (2), collecting terms of 0(dt), dividing by 
d^xdt, and taking the limit dt o, yields the bubble transport equation ; 
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8t|)/9t + V-V^i|j + a-v^ip + vd\p/dz = S-E^I|) 



(7) 



v^ere 



E. =: V ‘a. + 9v/9Jl 
tv 



( 8 ) 



and vAiere 



-» < x,y,z < “ 




< o>, and o < £ < “ 



A simpler (bat less rigorous) derivation of equation (7) is obtained by 
recxognizing that aij), and vip are the fluxes of bubbles in position space, 
velocity space, and radius space, respectively. Then 



is the net flux of bubbles into d^x by virtue of "streaming" in seven-dimensional 
space, (r,v,£) . Since the time rate of change of i|) in d^x is due to bubble 
streaming into d^x, as well as to other soiorces, S, we have 9iJ;/9t = F + S. Using 
the identity (ai|)) = and the fact that r and v are independent 

variables (V^*v = o) yields equation (7) . 

The term E^ij) in the bubble transport equation is generally nonzero due to the 
functional dependence of acceleration, a(r,v,£,t), and radius change rate, 
v(r,v, £,t), on the bxabble velocity and radius, respectively. This is uniasual since 
in transport equations normally encountered, the term corresponding to is 
generally zero (Chapnan & Cowling 1964, p.46 and p.322) . Mathematically, ^ o 
causes the volume element 



F = - { V^- (vi|)) + + 9(vi|))/9£} 



d^x' = [l+Z^dt+0(dt^) ]d^x 
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to change appreciably in the tintie interval dt. Physically, the. term on 

the right hand side of (7) acts as a psexado source, or sink, depending upon 
the functional form of a and v. 

For exaitple, suppose a is proportional to (-v) . Then V^*a < o, and the 

bubbles are decelerated into a smller volume, d^v' = (1+V^*^t)d^v. This 

increases the density so -(V^*a)i|) acts like a "soturce" in (7). Similarly, 

suppose V is proportional to v Jl so that dv/dZ is proportional to v . For 

z z 

V > o the bubble radii are increasing and becoming more spread out over a 
z 

larger interval, dJl' = (l+(3v/3£.)dt)dJl. This decreases the density, so -(3v/3£.)iJ^ 
acts like a "sink" in (7) . 

The second term on the right hand side of (7) , S(r,v,i,t) , represents 
distributed bubble sources and sinks. In addition to external sources and 
sinks, this term can include the effects of phenomena such as bxjbble scattering 
off small scale turbulent eddies, bubble-bubble interactions, etc. For 
instance, vdien the scale of turbulent eddies is the order of bubble radius, they 
will be referred to as "small scale," and bubbles entrained by the fluid may 
scatter off of them. To handle this phenomenon one could define a turbxilent 
scattering cross section, Z^, such that 

Zg (r ,Vj^^v, , t) v^ijj (r ,Vj^, 2,^ , t) d V^d2,^d‘^T 

represents the probable rate at which bihbles in d r about r at time t are 
scattered by turbulent eddies from d^v^^ about v^^ and d2j^ about 2,^^ into dV 
about V and d2 about 2, where v^ = Iv^^l is the sp^eed. Then 

S(r,v,2,t) = J 

(9) 

-v\jj(r,v,2,t) J" Z^(r,v-*Vj^fZ->Z^,t)d^v^dZ^ 
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represents the net rate at vMch bxjbbles are introduced into d^t. Specification 
of E requires detailed knowledge of the turbulent field as well as a model for 
bubble-eddy interaction. 

Substitution of equation (9) into (7) resiilts in an integro-differential 
bubble transport equation. If bubble-bijbble interactions are included, the 
b\±ble transport equation is nonlinear . as well (Chapman and Cowling 1964, p.63) . 
Below the surface in the upper ocean the average separation between bubbles is 
about 1000 times the bubble radii (Medwin 1970) , but very near or at the surface 
wave action can cause higher bubble densities. In addition, the energy contained 
small scale turbulence appears to increase with wave action but to decrease with 
depth (Shonting 1968) . lb avoid at this time obvious analytical corplications , 
the bubble density will be assumed low enough so that bubble-bubble interactions 
are negligible. It will also be assumed that the energy contained in small scale 
turbulence is negligibly snail, so the bubbles are partially entrained by the 
flxoid without abrupt changes in velocity or radius. These assumptions should be 
valid except very near the surface in heavy seas. 

Suppose the source, S(r,v,Jl,t), is independent of rp or, at most, is 
proportional to \p, 



S = Si + E^ip (10) 

This mathematical simplification yields a linear, first order, partial differ- 
ential bubble transport equation vbose formal solution is outlined in the next 
section. It is worthvbile to note that the simplifying assuitption (10) is not 
a severe limitation for many problems involving gas bubble transport in the 
upper ocean. All external bubble sources such as surface waves, rain and snow 
nucleation, organic decay, and photosynthesis can be included, as well as simple 
models of bubble scattering, absorption, and creation that have the form (10) . 
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III. SOLUTION BY METHOD OF CHARACTERISTICS 



The first order partial differential equation (7) is equivalent to the 
following set of eight simultaneous, first order, ordinary differential 
equations with initial conditions: 



dr/dt = V, -oo < x,y,z < «>, 


r(t^) 


= r = (x ,y ,z ) 

O 0^0 o 


(lla-c) 


dv/dt = a, -<» < V ,v ,v < «>, 
' X y z 


v(y 


= V = (v ,v ,v ) 
o xo yo zo 


(lld-f) 


dZ/dt = V, o < Jl < <», 


My 


o 

II 


dig) 


dip/dt = o s ^ 




0 

II 


(12) 



The equivalence between these "characteristic equations" and the bubble trans- 
port equation is most lucidly doronstrated with a geometric argument given by 
Garabedian (1964, p.l8) . The parametrically-represented cxirve C in nine- 
dimensional space, 

S(t) H [t,x(t) ,y(t) ,z(t) ,v (t) ,v (t) ,v (t) ,i!.(t) ,t);(t)] , -«> < t < <», 

X y ^ 

obtained by integrating the characteristic equations (11) and (12) , is called 
a "characteristic curve." Components of the initial point, S(t^) , are called 
"characteristics." The family of characteristic curves, 

{S(t) :-oo < t < 00 , S(y ep}, 

arameterizedby the set P of accessible characteristics, forms a hypersuTface 
in nine-dimensional space, tfi (t;r^,v^,£^t^) , that satisfies the bubble transport 
equation. 

The family of characteristic curves is essentially parameterized by the set 

(r ,v ,Z ) of accessible initial conditions since t and \p (r ,v ,Z ) , the 
' o o o o o o o o 

initial distribution, are generally fixed for a given problem. Since a, v, 
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and S are generally functions of r,v,Z, and t, the solutions to (11) and (12) 

are not only functions of the parameter t, but they also depend v^on all the 

initial conditions (r ,v , £ ) : 

o o' o 



r = r(t;r ,v ,t ) 

o' o' o' o 


= r(t) 


(13a) 




= v(t) 


(13b) 




= £(t) 


(13c) 


and = iHt;r^,v^,£^,t^) 


E ij;(t) 


(13d) 



\4iere the short notation on the right is used for brevity. 

The characteristic equations (11) are sinply the bubble dynamics equations. 
The bubble dynamics mcdel used to construct a and v (see ^4) , as well as the 
characteristics determine a "trajectory" (r (t) ,v(t) ,£ (t) ) vMch is 

the projection of a characteristic curve frcm nine-dimensional space onto 
seven-dimensional space. Bubble trajectories in (x,y,z) space are likewise 
projections of characteristic curves. 

A. Lagrangian and Eulerian Representations 

The parametric representation (13) , c±>tained by integrating the character- 
istic equations simultaneously, describes the continuum fran a Lagrangian frame 
of reference. In this description, bubbles are labeled by group according bo 
their characteristics (r^,v^,£^). This moving frame follows a specific group, 
initially in d^T = d^r d^v d£ about (r ,v ,£ ) at t , along a characteristic 
curve to d^x(t) about (r (t) ,v(t) ,£ (t) ) at t. At any time t, i|;(t)d^T(t) is the 
number of bubbles, belonging to the group labeled that are found 

in d^x(t) about (r (t) ,v(t) ,£ (t) ) . 
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On the other hand, the distribution ijj(r,v,£,t) , expressed vdth r,v,£, and 
t as independent variables, more conveniently describes the continuum frctn an 
Eulerian viewpoint. The Eulerian frame is fixed, and attention is focused on 
the particular volume d^x about the point (r,v,Z) . The instantaneoxxs bi±>ble 
density ij;(r,v, )l,t)d^T is considered without regard to the characteristics of 
bubbles in the volume. 

. Formal Solution 

If the source, S(r,v, Ji,t), has the form (10) , then (7) is analogous to the 
Boltzmann transport equation written for neutral particles (e.g., photons or 
neutrons) or charged particles (e.g., electrons) in a pxirely absorbing medium 
with distributed sources (Case and Zweifel 1967, p.31 and p.247). For exarrple, 
neutrons will stream along their characteristic curves in (x,y,z) space, which 
are straight lines, until they are absorbed or until they escape the medium. 
Similarly, the bubbles described by (7) stream along their characteristic 
curves in (r,v, £) space until they disaj^ar (Z o) or until they escape the 
medium. However, the characteristic curves defined by (11) are generally not 
straight lines in (r,v,2.) space or any of its subspaces, such as (x,y,z) . 

With S given by (10) , the formal solution to (7) is obtained by integrating 
( 12 ): 



+ rSj^(r(t') ,v(t')£(t') ,t')exp[Tj^(t,t')-T 2 (t,t')]dt' (14) 



vbere 



Ti(t,y = ^ Z^(r(t'),v(t'),Jl(t'),t')dt' 



and 



T 2 (t,y H ^ Z^(r(t'),v(t'),£(t'),t')dt'. 
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All integrals in this expressicai, and all similar integrals appearing elsewhere 
in this paper, are xmderstood to be line integrals along the characteristic 
curve that passes through (^q'^o'*'o^ (r,v,£). 

To obtain the nunber of bubbles in d^T about (r,v,£) at time t, (14) is 
multiplied by this volume element. However, a consistent representation for 
the right hand side of the equation requires use of the Jaccbian, 



Y(t,f) = 8(r(t),v(t),£(t))/9(r(t'),v(t'),£(t')), 

vdiich transforms the volume element from one point to another on the same 
characteristic curve: 

d^T = Y(t,t')d^T ' . 

Using the expression derived in the appendix, 

Y(t,t') = exp[T2 (t,t') ] , (15) 

the solution can be rewritten in the convenient form 

i(^(r,v,«.,t)d'^T = exp[Tj^(t,t^)]ij;(r^,v^,£^,t^)d‘^T^ 

pt (16) 

+ / exp[T^(t,t')]Sj^(r' ,v' ,£' ,t')d'^T'dt' , 

v^ere the abbreviations r = r(t;r ,v ,£ ,t ), r' = r(t';r ,v ,i ,t ) , etc., have 

o o' o' o o o o o 

been used. Fran this expression it is apparent that the nunber of bubbles in 
d^T about (r,v,£) at time t is a sum along the characteristic curve that passes 
through (r,v,S.) at t. It is the nunber originally in d about 
t^ plus (minus) those added (subtracted) by sources (sinks) along the character- 
istic curve between (?/V, £). If ^ o in (10) then the result 

is amplified (attenuated) by the factors exp[T^] because of the creation 
(absorption) rate along the path. For t = t^ + dt, (16) expanded in a Taylor 
series to 0(dt) yields equation (2), as it must. 
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Equations (7) and (16) are equivalent statements about the physical model 
outlined at the beginning of section tvo. They formally represent relation- 
ships that exist between the bubble distribution, single bxibble (^niamics, and 
the bubble sources. As such, these equations can be used for calculaticns to 
investigate relationships among the three in real situations. For exaitple, 
if ai^ two are specified, properties of the third can be inferred, either 
directly or indirectly. 

Prior to using (16) for calculations, a bubble dynamics model (a & v) must 
be specified SD that equations (11) can be integrated to obtain the ctiaracter- 
istic curves. A fair amount of information is available concerning single 
bubble dynamics, and in the next section some of the inportant points are 
summarized. 

The distributed sources and flxoid properties, as well as the appropriate 
bubble dynamics model, are dependent on the particular situation being analyzed. 
In section five the problem of gas bubble transport in the upper ocean is used 
as the basis for sate general remarks concerning bubble sources and fluid 
velocity fields. Then equations (7) and (11) are written for a simple ocean 
model, and seme calculations are performed to demonstrate how (16) can be used 
to evaluate the bubble distribution vAien the source function and bubble 
dynamics have been specific. 

It is significant to note that both (14) and (16) represent the distribution 
from a Lagrangian view point, i|) (t;r^,v^,£^,t^) . Transformation to the more con- 
venient Eulerian expression, \p(r,v,i!,,t) , requires the inverse of equations (13), 

namely r (t ;r,v,Jl,t), v (t ;r,v, £,t), and Z (t ;r,v,Jl,t). In general, 

O O O O o o 

equations (11) will have to be integrated numerically to establish a mapping 
between (^q'^o'^o'^o^ (r,v,2.,t), but often approximations can be made that 

will yield analytical expressions. 
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If S is a more corplicated function of tp than (10) , the formal solution (16) 
represents an integral equation for \jj. When analytical methods aimed at solving 
(7) or (16) fail, the Neumann series solution to the integral equation can be 
used to generate numerical results, if the series converges fast enough. 
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IV. SINGLE BUBBLE DYNAMICS 



In this section general expressions are derived for a and v, and some known 
models of gas diffusion and drag are sunmarized for gas bubbles in liquids. 

For the purpose of discussion in sections four and five, consider an inertial 
cartesian coordinate system with the z-axis vertical v^ward and the i origin at 
the surface of the liquid. Let v be the bubble velocity and V(r,t) be the 
transporting fluid's velocity relative to this frame of reference. 

A, Acceleration 

If a bubble of volume a, containing fluid of density p, were completely 
entrained (v = ^) by a fluid of density p^, it would experience the same 
force, p^adV/dt, as would transporting fluid enclosed in the same volume. For 
p 7 ^ p^, partial entrainment occvirs, and the bubble maintains a relative 
velocity, u = v-^, with respect to the transporting fluid. Hence, the bubble 
experiences a drag force, vdiich is discussed shortly in the context of 
gas bubbles in liquids . 

As the bubble moves relative to the transporting fluid it "drags" with it 
an amount of fluid having a volume equal to seme fraction, g, of the bubble 
volume. By Newton's third law, this adhering itass provides a supplemental 
reaction term, - ( g a p^ du/dt) , that acts to increase the b\Jbble ' s effective 
inertial mass. In addition, there is a net buoyant force, ag(p^-p)k, where 

A 

k is a unit vector vertically upward and g is the acceleration of gravity. 
Conbining all of these forces in Newton's second law yields the following 
repression for acceleration: 

a = ?j/(ap^(g+p/p^)) + (1+g) (dV/dt)/(g+p/p^) 

+ (l-p/Po)gk/(g+p/Po) (17) 
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in v±iich 



d^/dt = 3^/3 1 + (v*vp^ 



(18) 



B. Drag on Gas Bubbles in a Liquid 

Because of its theoretical and practical iitportance, the notion of gas 
bubbles in liquids (e.g., air bubbles in water) has been actively studied 
(levich 1962, ^80). Ihe regimes of bubble motion are classified according to 
Reynold's number, vdiere n is the liquid's viscosity. Cbservations 

of different-size bubbles rising in various liquids indicate that small 
(Rg < 1) and medium-size (1 < < 700) biibbles maintain a spherical shape, 

v^ereas large (R^ > 700) bubbles deform to flattened ellipsoids and very large 
bubbles (R^ > 4500) are unstable and tend to break up. 

For small bubbles (£• ^ 10 Oy = 10 “‘'to in water) the drag is viscous in nature 
and is given by 



with K = 47 t for a perfectly clean bubble having a mobile two fluid interface 
at its surface. However, most liquids such as water contain "surface-active" 
materials that coat the bubble and destroy the mobility of this interface. In 
this case the relative velocity of the transporting fluid goes to zero at the 
bubble's surface, and the bubble behaves like a solid sphere vdiere drag is 
given by (19) with k = 6tt (Levich 1962, 70 and 81) . 

The flow past medium-size bubbles (lOOy < I < 2000y in water) is separated, 
with the separation region occupying an area s^^ on the downstream portion of 
the bubble's surface. Up to the separation point resistance that acts on the 
bubble is viscous in nature with a contribution to the total drag given by (19) 




(19) 
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with K = 12 it. Past the separation point and into the bubble's wake the flow 
is characterized by turbulent motion with a contribution to the total drag 
given by 

= -O.SK^p^Sj^uu, (20) 

vdiere is a drag coefficient with 0.65 > > 0.55 for 200 < 1000 

(Levich 1962, 80 and 82) . 

In the absence of surface active naterials, is very small 

(Sj^ ^ ^ flow past the bubble is essentially unseparated 

and the drag is given by (19) with k = 12ti. On the other hand, when an area 

s^ on the bul±>le is covered by a monolayer of surface-active material, the 

relative fluid velocity in this region is zero, and flew separation occurs 

there. In this case viscous drag (19) is accorpanied by the form drag (20) 

with s, = s , and the latter dcaninates the fomer vhen 
1 o 

s^/(47t)i2) ; 28 (n/p^)V{g)l^) 

(Levich 1962, ^ 82) . For exaitple, \dien i = 500y, only 2% of the bubble's 
surface need by covered by a monolayer of this material before the drag force 
becanes quadratic according to (20) . 

C. Radius Change Eate 

For gas bubbles moving in a liquid such as water, radius changes are caused 
primarily by changes in pressure, due to changes in depth or surface tension, 
or 1^ gsis diffusion across the bubble surface. Spherical bubbles (R^ < 700) 
that contain n moles of ideal gas at tertperature T and pressure P are described 
by 

4ttJiV 3 = nRT/P, (21) 
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where R is the gas constant. For isothermal processes, (21) is readily 
differentiated to yield 



,, = = 1 / dn/dt _ dP/dt \ 

dt 3 \ n P /* 

Gas diffusion is included in the dn/dt term, v/hile corrpression 
and surface tension are included in dP/dt. Assuming quasi-equilibrium, gas 
pressure is the sum of atmospheric pressure, P^, srurface tension pressure, 2c/£, 
and fluid weight, vtiere -z ^ o is the depth and c, is the surface tension? 

P = P^ — p^gz + 2^/Z, (23) 

and 

dP/dt = -p gv - (2?/£2)d£/dt. (24) 

o z 



D. Gas Diffusion 

As results of the next section will indicate, gas diffusion is one of the 
more inportant factors affecting the bubble distribution. It is a caiplex 
phencmenon depending rpon many influences that include the type of gas in 

the b\±ble, the gas diffusivity, D, in the transporting liquid, the gas con- 
centration, C, in liquid contacting the bubble's surface, the gas concentration, 
C^, in liquid far away frctn the bubble, the presence of surface-active materials 
in the liquid and on the surface of the bubble, the flew field around the bubble 
and the bijbble radius. To sirtplify the discussion here it will be assumed that 
gas inside the bubble is cotrposed of one type of "average" molecule having a 
single diffusivity, D. 

If the Peclet number is large, P^ = (u£/I)) > > 1, the gas concentration 
gradient is confined to a thin boundary layer on the surface of the bubble, and 

diffusion between C and C occurs across a small distance d 5,/P Hence 

00 © 

(Levich 1962, H 14 and 72) 
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dn/dt = -G«,^ (C-C ) , 

oo 



( 25 ) 



v^ere the cxsefficient G is generally a function of %, u, and D. For specific 
^plications in vhich individual corponents of gas must be monitored, equation 
(25) will have to be generalized. 

typically (n/p)/D ''^10^, so there is a wide range of bubble radii for vhich 
the P^clet number is large (P^ > > 1) but the Reynold's number is small (R^ < 1) . 
For that case (small bubbles) , 

G = 8[(tt/ 6) (DuA)]l/2 (26) 

in the absence of surface-active material (Levich 1962, 72 and 91) . In the 

presQice of surface-active material, the flow field around a small bubble is 
like that around a solid sphere, and G becomes (Levich 1962, ^ 14) 

G = 8(d2u/«,2) 1/3 . (27) 

A representative value for u is the terminal speed, u,p, of small bubbles rising 
in a quiescent fluid. Equating the drag (19) to the buoyant force, 4TrJl2gp^/3, 
vhich assimnes p/p^ < < 1/ results in the expression 

u,j, = (4 tt/3k) (gpQ«.Vn) , (28) 

vhere k = 4tt in the absence of surface-active materials, and k = 6tt in their 
presence. 
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At moderate Reynold's ntnribers (mediijin-size bubbles) the flow is separated, 
but the region of separation, s^^, is generally small. If gas diffusion across 
s^ can be neglected, (25) represents the gas diffusion with G equal to (26) 
multiplied by (Levich 1962, ^ 91) . 

If liquid in contact with the bubble surface is assumed to always be 
saturated with gas at the pressure P (in atmospheres) prevailing inside the 
b\±fole, then C = KP, where K is the absorption coefficient (kg*m“^’atm“^) for 
the liquid-gas system. Defining the partial pressure f = C (in atm) of gas 
in the transporting fluid, (25) can be ran/ritten 

dn/dt = -GK«.2(p-f). (29) 

Under laboratory conditions, Wyman, et al (1952) , measured the rate of 
solution and the change in conposition of air bubbles in stirred sea water, as 
a function of depth. These experimental results follow very closely 

dn/dt = -4Tr6Jl2(p-f) (30) 

with 6= 4.4 X 10“^ moles m“^sec“^atm“^ . Using equations (26) and (27) with pub- 
lished values of D and K for ojQ^gen and nitrogen in water, estiimtes of GK are 
found catparable to 4 tt 6, but generally larger. The presence of surface-active 
materials that retard gas diffusion and reduce the effective absorption 
coefficient across the bubble's surface may have contributed to the discrepancy 
(Fox & Herzfeld 1954) . 
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V. APPLICATION TO THE UPPER OCEAN 



The purpose of this section is to briefly discuss application of bubble 
transport theory to the description of gas bubble populations in the upper 
ocean. Because it involves a variety of interrelated phenotnena pertaining to 
different disciplines, this is a difficult subject to analyze. Also, there 
is only limited information available upon vMch to base mathematical models 
of bubble sources, fluid velocity fields, gas diffusion rates, etc. However, 
it is the author's opinion that bubble transport theory is a useful tool for 
interpreting data and relating phenonena that affect the distribution. 

To begin, it is appropriate to summarize sane observations concerning 
bubble sources, S, and fluid velocity fields, in the upper ocean and to 
iteke a few remarks concerning analysis. Ihen, to demonstrate an application 
of bubble transport theory, some sample calculations are performed and 
interpreted. 

A. Bubble Sources and Fluid Velocity Fields 

Physically one might consider bubble sources to be segregated into three 
categories according to depth: Those on the ocean floor, those concentrated 

at the surface, and those distributed throughout the medium, principally near 
the sxirface. Mathematically these classifications are convenient for analysis 
since sources in the first two categories can be incorporated into (7) and (16) 
either as boundary conditions or as plane sources. 

The first category includes sediment-initiated bubbles vfeich often contain 
organic gas. For exantple, McCartney and Bary (1965) measured relatively large 
bubbles with radii between 450y and 800y caning off the botton of Saanich Inlet, 
British Columbia, and speculated that these bubbles contained methane. 
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In the second category are bubbles generated at or very near the surface 
by wind and wave action, precipitation, and captured aerosols. The subsurface 
measurements of Blanchard and Wbodcock (1957) , Glotov, et al (1962) , and 
Medwin (1970) imply that most of the bubbles associated with breaking waves 
have radii less than 200u. On the other hand, experiments conducted by Monahan 
and Zietlcw (1969) , in which a small volume of seawater was poured into a tank 
of the same, indicate that the cloud of bubbles formed by the descending plime 
of mixed air and water contains an appreciable nuirber of 200p to lOOOp bi±)bles 
as well. In addition to wave-generated bubbles, rain and snow nucleate bubbles 
at the surface vidch have radii predoninantly less than 50p (Blanchard and 
Woodcock 1957) . To explain the excessive sound attenuation and scatter 
observed in near-surface water at high frequencies (f > 60kHz) , Medwin (1970) 
speculated that continental aerosols may be the primary source of small bubbles 
(Z < 50yu) near the surface. 

In the third category are bubbles initiated within the volume of the upper 
ocean. Internal wave action, radiation, and biological activity, such as photo 
synthesis and marine life gas emissions, have been postulated to be responsible 
for many of the smallest bubbles (Z < 60p) . Also, seme of the volume sources 
may be associated with pressure fluctuations caused by wind-driven surface 
waves. Pressure changes and turbulence may nucleate small bubbles at distri- 
buted cavitation sites, vhile the vertical fluid velocity and turbulence 
created by breaking waves introduces larger bubbles belcw the surface. 

When a wave breaks and air is rapidly mixed with water, the bubbles fonred 
are entrained by a strong vertical velocity field; these bubbles have been 
observed as far as three wave heights below the surface (Kanwisher 1963) . As 
the strong downward currents decay, the bubbles tend to rise under buoyancy 
with motion modified by the remaining fluid velocity field (ref. equation (17)) 
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The largest bubbles rise most quickly and reach the surface before going into 
solution v^ile the smallest bubbles succurb to gas diffusion (LeBlond 1969) . 

As a result/ the most prevalent bubble within a few meters of the surface should 
have a radius between the two extremes, an assertion supported by the data of 
Medwin (1970) and Glotov, et ^ (1962) \dio measure distribution maxima at 90y 
and 60y, respectively. 

Shonting (1968) has measured the autospectra of particle motions in the 
upper ocean and shown that the dominant peaks are associated with the fre- 
quencies of the surface waves. The spectral energy decays exponentially with 
depth in a way that attenmtes high frequencies the most rapidly. 

Example 

Consider a one-dimensional, steac^ state ocean in vbich all velocities are 
vertical and in vbich all functions depend spatially on depth only. In this 
case the bubble transport equation reduces to 

v^ere^ (z,v,£) is the distribution function, S(z,v,£) is the distributed soxarce, 

V = V is the vertical velocity, a = a(z,v, £) is the vertical acceleration, 
z z 

v(z,v,£) is the radius change rate, and = 9a/9v + 9v/3£. 

The characteristic equations for (31) can be written with z ^ o chosen as 
the independent parameter. This choice is equivalent to dividing equations 
(11) and (12) by (11c) : 

dv/dz = a(z,v,£)A^, -«> < v < » & v(z^) = v^; (32a) 

d£/dz = v(z,v, £)/v, o ^ £ < “ & (Zq) == (32b) 

d(j)/dz + ('Z^/v)({) = S(z,v,£)A', o ^ <}i < “ & 'I>(Zq) = 
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Sijnultaneous integration of equations (32) yields the family of characteristic 



curves 



{z, v(z;z^,v^,J!,^) , £(z;z^,v^,Jl^) ? z ^ o} 



(34) 



and bubble trajectories in (x,y,z) space are vertical straight lines. 

A formal solution analogous to (14) is readily obtcdned by integrating (33) 
along characteristic curves (34) . The form analogous to (16) , 

V(J) (z,v,£)dvd£ = v^<}>(z^,v^,£^)dvjd£^ + / S(z',v',£')dz'dv'd£', (35) 

is derived using the Jacobian, 

= a ' iv'^l,*) ^ ^ ^ ^ ' (36) 

' Jz» 

Expression (36) follovs from the differential equation 



dy/dz = (E^ - (a/v2))y 



(37) 



obtained when characteristic equations (32) are catibined in a manner ..similar to 
the general scheme detailed in the appendix. Of course, all integrals 
appearing in this section are taken along characteristic curves (34) , vhere 
abbreviated notation v' e v(z';z^,v^,£^) £" 5 £(z”;z^,v^,£^) , etc. has been 

used in (35) and (36) . 

Physically equation (35) states that the flux of bubbles, V((), at the point 

(z,v,£) , equals the flux of bubbles, ^.t the point ( 2 q>v^,£^), plus (minus) 

the bubbles added (removed) by distributed sources (sinks) along characteristic 

curve (34) . Further physical insight can be gained by reflecting on a situation 

where there are no distributed sources, S = o, so that <}>(z;z ,v £ )dvd£ = 

^ o o o 

(v^/v) (J) (z^,v^,£^)dv^d£^. If there is no fluid velocity, bubbles rise under 
buoyancy at terminal velocity and may shrink and slow down or expand and speed 
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up (LeBlond 1969) . Those belonging to the Lagrangian group v^iose 

velocity and radii increase vath z have a density that decreases as z increases 
because the bvibbles at z + dz are running away from those at z. On the other 
hand, those v^ose velocity and radii decrease as z increases have a density 
increasing with z since the bubbles at z are piling up on those at z + dz. 

1. Characteristic Curves 

A first step toward using (35) to calculate the bubble distribution is the 
specification of a(z,v,£) and v(z,v,£) and the ccrrputation of characteristic 
curves (34) . A usefiil expression for the radius change rate is obtained by 
substituting (21) , (23) , (24) , and (29) into (22) and solving for d2./dt. Then 
(32b) becomes 

dZ/dz = v/v = { (V3)-(GK/4 tt) (KT) [D(l-f)-z + (r/Jl) ]/v}/[D-z + (2/3) (r/Z)] (38) 

vbere (-z) is the depth, f is the partial pressure (in atmosEberes) of gas in 
the water, 

D 5 p ^/p^g = 10 meters 

is the depth of sea water equivalent to one atmosphere, and 

r = 2c/p^g = 1.47x10“^ maters^ 

is a surface tension coefficient in water. Note that the r/Z terms are 
negligible unless the bubbles are very small (Z < 30y) . 
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Acx:eleration is given by (17) with the drag force given by (19) or (20) . 

If the fluid velocity is negligible (V = o) , the bubbles rise very nearly at 
their terminal velocity, v^. This is seen in Figure 1 vhich was obtained by 
simultaneously integrating equations (32) ramerically* for bubbles of various 
initial radii, starting from rest at = -20 meters. Small bubbles 
{I < lOOy) , and medium-size bubbles (lOOy < I < 2000y) with viscous drag domi- 
nant, rise at terminal velocity v^ given approximately equating (19) to 
the buoyant force; 

(39) 

where k' = 1 and k' = 2, respectively, and vhere 

a 5 (9/2) (ti/Pq) = 4.37x10”^ meters Vsec 

in water. When medium-size bubbles have a sufficient fraction of their surface, 

e = s^/(4ttJ12) , 

covered by a monolayer of surface active rraterial, turbulent drag dominates 
viscous drag because flow separation occurs at the boundary of s^. They rise 
at terminal velocity v ^2 given approximately by equating (20) to the bixayant 
force; 

v^2 = (g/5) (40) 
vhere K = (3/2) e = 0.9e. 



*A11 numerical integrations alluded to in this section were carried out xasing 
gaussian quadrature or optimal fourth rank Runge-Kutta (Ceschino and Kuntzmann 
1966, p. 67) . 
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Datta, et al, (1950) has averaged the results of several different observers 
to produce an eirpirical curve for the terminal velocity of air bubbles in v/ater 
as a function of their radii. Values produced by (40) agree with this 
empirical data for 300y < £ < lOOOy if the surface fraction, e, is reduced mono- 
tonically as the radius, %, increases. This leads one to speculate that as the 
radius and terminal velocity increase, rapid liquid motion at the siirface either 
washes away surface active itolecules or carries them to the rear \diere they form 
a ccmpact, saturated monolayer (Levich 1962, p.447) . 

In general the terminal velocil^ of a rising bubble is not given simply by 
(39) or (40) during its entire lifetime. Many bubbles shrink or ejpand as they 
rise (LeBlond 1969) and thereby transition from the regime of predominantly 
turbulent drag to the regime of predominantly viscous drag, or vice-versa. The 
treatment of separated flew in Levich (1962, 82) suggests that the sum of (19) 

and (20) should be a useful model for the drag on medium-size bubbles in the 
transition region. Using this model in (17) , (32a) can be written 

dv/dz = -[ (5/6)/£] [ (v-v_) (v+v^)/v] , (41) 

where 

v± H + (1+4q2)1/2]/(2Q) (42) 

with 

and 

I = [(K'a)2/(ga]l/3 

Numerical integration of (38) and (41) demonstrates that the bubbles released 
from rest rise with terminal velocity v^ = v_ (figure 1) . This approach appears 
to have some merit since (42) for v_ yields values close to the empirical curve 
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by Datta, et al (1950) and since the model shews the proper limiting behavior; 
Pbr SL < SL, viscous drag dominates turbiilent drag in (41) and 



v_ = v^j^ - v^^{q2 + 0 ( 0 "+) } 



for 1 > i, turbulent drag deminates in (41) and 



v_ = v ^2 + 0 (Q“ 2 ) }. 



Conputation of characteristic curves (34) in a quiescent fluid (V = o) is 
simplified because the bubble velocity is always very nearly the terminal 
velocity, vsv^(z;z^,il^) , and is essentially independent of v^. Either (39), 

(40) , (42) , or aipirical values for v^ can be used in (38) to obtain i{z}z^,Z^) . 
If (39) or (40) is used and surface tension is neglected, closed analytical 
expressions result: For v = v^^^, r = o, and f constant. 




D-z \ 
o) 




(43a) 



with an inverse transformation 




(43b) 



for V = r = o and f constant. 





(44a) 
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with an inverse transfontation 



^•ip2 (ZqJ z / £) . (44b) 

Otherwise, (38) must be integrated numerically. 

Tb evaluate accuracy of approximate analytical results (43) and (44) , 
and are plotted vs depth in figures 2,3, and 4 and compared to , 

obtained by numerically integrating (38) with v = v_ given by (42) . For the 
sake of exartple, f = 1 and KrGK/4ir = 10® meters/sec (Wyman, et al 1952) are 
used in all three models, and il^^(z), ^rp(z) are represented for 

various initial radii, , at z = -20 meters. Note that is a good 
approximation to for 30p < i < lOOp (viscous regime) and that £^2 closely 
approximates for £ > 400y (turbulent regime) . Neither (43) nor (44) is 
valid for £ < 30p because they neglect surface tension, and neither is accurate 
in the transition region lOOy < £ < 400y. 

An exaitple of the phenonenon analyzed by LeBlond is evident in Figures 
2-4. As a b\±)ble rises under buoyancy it tends to expand as the hydrostatic 
pressxire decreases, but it also tends to contract as gas diffuses across the 
bubble's surface. If a bubble's radius is less than some critical radius, 
it rises too slowly, shrinks monotonically, and disappears before reaching the 
surface because gas diffusion dcminates hydrostatic expansion. If its radius 
is greater than another critical radius, £ , it will expand monotonically as 
it rises because hydrostatic ejq>ansion dominates gas diffusion. For an initial 
radius between £^ and the bubble's behavior is not necessarily monotonic; 
it may or may not reach the surface. 

Bubbles with radii less than £j^ will rise a finite distance before 
iisappearing, and carputed values of £^ and £j^ depend upon the gas diffusion 
todel used as well as parameters like z^. Examination of £,p in Figures 2—4 
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indicates that 400y and 'v 250y, while a lOOy bubble rises about 1.1 
meters. Unfortunately there is usually insufficient Infonnation available 
for an accurate specification of the gas diffusion phenonenon. For exartple, 
in (29) the gas diffusion rate depends not only upon the partial pressure f (z) , 
but also upon the parameters G and K viiich are determined by the type of gas 
and the presence of surface active materials. This data is often 
unavailable for the ocean. Since characteristic curves and bubble distributions 
(ref. Figures 5 -13 are sensitive to the gas diffusion rate, specification of 
the unknown parameters is inportant, and it may require a sequence of educated 
guesses refined fron ccrrparisons of cosnrputed and observed distributions. 

If the fluid velocity ^ is not negligible and the problem is multi- 
dimensional, the conputation of characteristic curves is more difficult. In 

general, equations for x,y,v ,v , and v must be integrated simultaneously 

X y z 

with (38) . However, if d^/dt given by (18) is small caipared to the other terms 
in (17) , a useful approximation (that eliminates the three characteristic 
equations for velocity) might be 



v = V + v^. 



where v^ is the terminal velocity of the rising bubble. 

2 . Bubble Distribution 

Since bubbles rise in a quiescent fluid at very nearly their terminal 
velocity, 

(|)(2,v,il) = 6(v-v^) $(z,£) (45) 
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should be a good approximation vAiere is given either enpirically or by (39) , 
(40) or (42) . In this equation 



$(z,^.) = / 4>(z,v,)l)dv 



(46) 



is the "radius density/' the number of bubbles per unit volums per unit radius. 
Equation (35) can be integrated over v to yield 

§(z,A) = (v^/v^) §(z^,A^) OJiyajl) + (1/v^) / s(z',)i') 0)i'/3S,)dz' (47) 



with 



s(z',il') = S(z' ,v' ,Jl')dv' , 



(48) 



v4iere d)l = (32, /3)l)d)l and dn' = (32.'/32)d2, have been used. The bubble density, 



$(z) E / |(z,Jl)dJl, 
'o 



(49) 



which is the total number of bubbles per unit volume at depth z, regardless of 
radius, is given by (47) integrated over all radii. 



oo pz 



i(z) = (v^Vj,)§(z^,yd)i^ + 



d/v^)s(z',2,') (3£'/3il^)dz'dJl^ 



o J z 



(50) 



It must be remaribered that v^ and v^, as well as 2 and 2' , are functions of 2^ 
and that integrations in z are taken along the characteristic curves 



{z',v^(z';z^,2^) ,2^(z';z^,2^) :z^ _< z*^ z}. 



(51) 



where 2^ is obtained by integrating (38) with v = v^. 
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For most applications the radius density provides adequate information. 

To calculate this function for a given z and I, one needs di’/dl, and 
s(z',£') at points on the characteristic curve that passes through (z,Z) so 
the integral in (47) can be evaluated. Using notation from the appendix, 

= l/Y 77 (ZfZ')f (52) 

vdiere Y 77 satisfies the differential equation c±>tained by setting v = v^ in 
(38) and taking the partial derivative with respect to 

dy 77 /dz ^77 ”^77 (53) 

with the boundary condition 

Y77(z,z) = 1 

vdiere 

F- 7 (z;z^,£^) = {(2rA2) [£-(36RT/v^) (D(l-f)-z+r/£)] [3(D-z) + 2T/Z]~^ + 

[1+ (36RT/v^) [r/£2+(D(l-f)-z+r/£) Ov^3£)A.j,]]>[3(D-z) +2r/£]"l (54) 

in which (GK/4ir) = constant e 6 (V^man, et al, 1952) has been used. When the 
approximation v^ = v_ (42) is used, 

3v^3£ = (3v^^A) [1 + 4Q2]-1/2 - v_/£. (55) 



Ihe solution to (53) is 



Y 77 (z,z') = exp[7(z,z')] 



with 



J(z,z') E / F^^(z";z^,£^)dz", 



(56) 
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and it is obvious that 



Y^^(z,z') = Y7-7 (z,z^)/y^^(z',z^) . (57) 

Using (56) and (57) , (47) can be written 

§{z,SL) = [exp[-7(z,z^)]/v^] {v^§(z^,y + 

Pz (58) 

+ / s(z' ,Jl')exp[^(z' ,z^)]dz' } 

Once the source s(z,il) has been specified, (56) can be used to evaluate the 
radius density at points along any characteristic curve (51) obtained by 
integrating (38) . In the following discussion of saitple calculations, 

J(z,Jl) will refer to the radius densil^ when the approxiination = v__ (42) 
is used. 

If either approximation (43) or (44) is appropriate, then either 

= [(D-z)/(D-z')] (43c) 

or 

aJlV9^ = [(D-z)/(D-z')]l/2(ji/ji.]l/2^ (44c) 

is used in (47) to obtain the approximate radius density §^(z,l) or , 

respectively, at points along the characteristic curve (51) with 5,^ = or 
respectively. If gas diffusion can also be neglected (GK = o) then 
both approximations reduce to 

I = l^^{z}z^,l^) = I(D-z^)/(D-z)] (59) 

which can be used in (47) and (51) to obtain the approximate radius density 
(z , Jl) . 
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3, Numerical Results 



For the sake of exanple, consider the situation investigated McCartney 
and Bary (1965) who acoustically measured rather large gas bubbles ascending 
from the bottcm of Saanich Inlet, B. C., and inferred the botton source repre- 
sented in Figure 5 from the ascent velocities. If we model that ocean as 1-D 
and quiescent, then equation (47) can be used to infer the radius density 
throughout the volxjme. Since vol'ume sources produce mainly smaller bubbles 
that their measurements do not include, only the plane source f (z^,£^) 

(Figure 5 ) at = -197 meters is considered in this sample calculation, and 
(47) reduces to 

i(z,£) = (v^/v^)|(z^,£^) Onysj!.) . (60) 

Using f = 1, approximations ^3 ^ radiiis were coiputed 

at various depths between z = -197m and z = 0. f (z,£) (solid line) and 
§ 2 (z,£) (dashed line) are plotted in Figures 6-10 at the depths z =“172m, ~147r 
-122m, -97m, -47m, and Om. Although actual calculations do not yield 
precisely rectangular distributions for z ^ z^, $ and drawn 

as perfect histograms because the error introduced is negligible. 

While ^2 neglects gas diffusion and surface tension, ^ includes these 
effects. For the gas diffusion model used, all bubble radii for which 
^(Zo,£o^ is nonzero lie below the critical radius so these bubbles will 
shrink monotonically cis they rise and will disappear before reaching the surface. 
The distribution ^ is seen to shift toward smaller radii eis depth decreases, 
until even the largest bubbles have disappeared by the time z v -90m, so only 
§2 is shown in Figure 10* Since = o in the model used to calculate these 
bubbles necessarily expand as they rise, and the distribution shifts monotonically 
to larger radii. 
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Ihe bubble density is sinply the area under the histogram. As z increases 
frcm -197m the density of 4>(z,i!,) initially increases since flux is conserved 
and the bubbles are slowing down and piling up. However, nearer the surface, 
most smaller bubbles have succurribed to gas diffusion and the density starts to 
decrease toward zero. For ^^(z,^) the density is seen to decrease as z 
increases since the bubbles are accelerating and becoming spread out. 

The distributions $2 and are not shown because neither follows $ closely. 

$2 is the best approximation, but these bubbles do not quite shrink fast enough. 
Sane of the bubbles in contract monotonically v^le others contract and then 
expand. For the model used most of the bubbles originate at the botton with 
radii in the transition region, ^ ^ v/here neither approximations (43) 

nor (44) are valid. 

Interestingly enough, the experimental results of McCartney and Bary (1965) 
indicate the bubbles fron the botton of Saanich Inlet expiand as they rise 
according to the distribution $ 2 , vdiich ignores gas diffusion. Thus, it appears 
that diffusion is sonehcw inhibited across the surfaces of these bubbles which 
are catposed of "gas released frcm the highly organic, anaerobic sediment of the 
bottom, " possibly methane . Either the water was saturated with respect to the 
gas(es) in the bubbles, or these bubbles were coated with an active material that 
severely impeded gas diffusion. While it is significant to note that $ is cal- 
culated using GK/4-n = 6 (eqn. (30)) , v^ch is based upon experiments with 
relatively clean air bubbles (V^an, et al 1952) , and that f = 1 at all depths 
is not physically realistic, certain qualitative conclusions, based on a com- 
parison of $ and $2 observed distribution, are still valid. For example, 

gas diffusion appears to be negligible in this situation; and gas diffusion 
would profoundly affect the bubble distribution. The latter observation emphasizes 
the fact that results are strongly dependent on the bubble dynamics and ocean 
models used. 
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As the basis for a second sairple calculation, consider the supposition 
that the observed bubble distribution indicates the functional dependence of 
the volume source. To test this idea, consider a hypothetical volume source 



suggested by the measurements of ShuUcin (1968) and Medwin (1970) 



s(z,£) = e^^(£) 



(61) 



viiere 



L(£) = ( 



(£//) ^^*^for £ ^ £' 
0,for 0 ^ £ < £' 



(62) 



with h = 10 meters and l' = 20 u. Physically, a source decaying exponentially 

with depth might arise from photosynthesis, light sensitive bacteria, wave 

action, etc. ^fedwin (1970) has observed an exponential behavior for small 

bubbles (£ < 60y) , and both Medwin (1970) and Shulkin have observed that the 

near surface bubble density appears to fall off like some power of the bubble 

—3 5 

radius. The radius spectrum L(£)'^'£ ’ is chosen as an exairple. 

Assuming that the source is zero belcw z = -30 meters, and that f = 1, 
the radius density was computed using v^ = v^ in (38) and (58) . Figures 11 & 12 
sents log-log plots of f(z,£) (solid line) and s(z,£) (dashed line) vs radius 
for four depths, z = -29m, -20m, -10m, and Qm. A comparison of § and s indi- 
cates that the distribution follows approximately the functional form of the 
source for £ > lOOvi: 



h' = 6.5m. The distribution at z = -29 meters does not fit this pattern 
because it is too near the region assumed to be devoid of bubbles. 

For £ < lOOy the functional form of | is more complicated and does not 
necessarily follow (63) . Except at the surface, the curves are all concave 
down vhich indicates a deficiency of small bubbles. Gas diffusion and surface 



J(z,£) -V e' 



zA’£-3.5 



(63) 
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tension cause small biobbles to disappear much more rapidly than large b\±ibles. 
While the snail bubbles observed at a given depth originate primarily from 
local sources, itany of the larger bubbles observed originate at sources distri- 
buted well below that depth. However, for the model used, gas diffusion 
decreases rapidly with depth, so there is less discriminatd.on against the small 
bubbles in the near surface region. From Figures 11 and 12, 

z (64) 

These resiilts indicate that if gas diffusion is significant, the small 
bubble sources will probably fall off to a hi^er power in £ than the observed 
distribution, but the larger bvibble sources may have roughly the same 
£-d^>endence as the observed distribution. Furthermore, the depth depsendence 
of the bubble sources is likely to be fxjnctionally similar to that of the 
bubble popjulation. Of course these qualitative relationships between bi±ble 
population and bubble source are not likely to be applicable unless the model 

1 

I is valid. For example, any of the following can have a profound influence on 

I the bubble distribution and its relationship to the source: A significant 
fluid velocity, ^ o, a partial pressure, f (z) , that varies with depth, or 
surface active materials that impede gas diffusion (e.g., lower K in 
equation (24) ) . 

The distributions and ^2 were also cctrputed for this exairple but are 
not displayed. As expected, approximates ^ fairly well for 30y < £ < lOOy, 
and $2 is a good approximation to f for £ > 350y. 
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VI. CONCLUSIONS 



The main purpose of this paper has been to develop the general framework of 
bubble transport theory and then to demonstrate its application. The bubble 
transport equation and its solution relate the enserrble average behavior of 
a bubble population to dynamics of a single bubble and to properties of the 
transporting fluid. Among the most inportant aspects of gas bubble dynamics 
in liquids are bubble gas diffusion and drag, vjhich depend upon many different 
j^iysical parameters. For the upper ocean, information concerning these para- 
meters, as well as information about bubble sources and fluid velocity fields, 
is usually inccmplete. In this case, bubble transport theory cannot be 
used to predict bubble distributions; its primary utility is as an analytical 
tool for investigating relationships among the various carplex factors that 
affect the distribution and for refining models of the upper ocean. Future 
work will be concemed with analysis of more conplex ocean models. 
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APPENDIX 



It is shown here that the Jacobian, 



Y(t,t') 



3(r,v,&) 



(A-1) 



used to transform the volume element d^x between two points on a characteristic 
curve, 

d^rd^A = Y(t,t')d2r'dV'd£' , (A-2) 



satisfies the differential equation 

dy/dt = . 



(A-3) 



Using the boundary condition 
obtain 

Y(t,t') 



Y(t;t') = 1, 




this equation can be integrated to 




Z )dt" 
o 



1 



(A-4) 



where the indicated integration is along the characteristic curve passing through 
(rjvji,') and (r,v,Ji). 

Referring to the notation (5) , define the elements of a 7 x 7 matrix 

Y^j = 9qj^/9qj ' , i,j = 1,2,..., 7 (A~5) 

then the Jacobian is the determinant of this matrix, 

Y = det((Y^.)) = (A-6) 

(or j) 

vhere is the "minor" of the element Yj^j (Cj^j is y with the i rcM and 
th. 

column deleted) . 
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If the following elanents are defined. 






(A-7) 



R 5 = a , Rg = a^, and E 



= V 



then the characteristic equations ( 11 ) can be written 



dq^/dt — (q^^ ,q 2 f • • • /q>^ , t) i 1,2, •••,7 



(A- 8 ) 



Taking the partial derivatives of these equations with respect to the q j ' and 
using the cha i n rule, the set of differential equations for the y^j is obtained; 



dy.. 9R. 

J-3 = _i 

3t 3q 



^ ^ i ^ i^j> 1,2,. ..,7 



(A-9) 



’3qi, 



These equations are used to expand the derivative of (A- 6 ) taken with respect 
to the parameter t; 



dy 

dt 




(- 1 ) 



i+j dy. . 



3R. 



iic V-i 

dt ij Z— ^ j / . 3qj, ’kj 

i,j k 



y. . (A-10) 



Using (A- 6 ) , this expression can be rearranged to yield the following: 



dy 

at 



C — > i+j f 3R. 3R. -j 

“ ^ij 1 3^ ^ij ^ ^kj j 

if j kf^i 



(A-11) 



^3R. 



r i 3R. J 

3^2-^^"^^ ^kj ^ij 



^ i,k?^i 



k 

J 
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The term 






Y, . C. . 
k: ID 



easily be seen by expanding 



is zero since l^i. This can most 
each determinant C. . by minors of the 



Since 




-> 

a + 



H z 

t 



the result (A-2) is obtained. 



elements 



(A-12) 
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LIST OF FIGURE CAPTICNS 



Figure 1: 

Ihe velocity, v, of bubbles starting at rest and rising from an initial 
depth = -20 meters. The initial radixis is (in microns) , the terminal 

velocity is v^, and (z-z^) is the height above the release point. 

Figures 2, 3 , and 4 ; 

Bubble radius in microns as a function of depth in meters for various 
initial radii, (at z^ = -20 meters) , calculated using three different expres- 
sions for the terminal velocity, v^. Ihe solid line is Z^, obtained by Inte- 
grating ( 38 ) with V = v_ (equation 42), the broken line is Jl|j,^Cequation 43a), and 
the dashed line is Z^^ (equation 44a) . The gas diffusion model used (V7yman, et al 
1952 ) assumes a uniform partial pressure of one atmosphere, f = 1, and 
RTGK/4 t 7 = 10"*^ meters/second. 

Figure 5 ; 

The distribution, o(z^,£,^), in relative units, of bubbles rising off the 
bottom of Saanich Inlet, B. C. (IfcCartney and Bary 1965) where z^ = -197 meters 
and the bubble radius, Z^, is in microns. 

Figures 6,7,8, and 9: 

The distributions $(z,£) (solid line) and f^(z,Z) (dashed line) in relative 

units at the depths -172m, -147m, -122m, and -97m where bubble radius, 

Ji, is in microns. $ neglects gas diffusion while $ assumes a uniform 

“6 

partial pressure of one atmosphere, f = 1, and RTGK/4 tt = 10 meters/ 

sec (Wyman, et al 1952) . 

Figure 10 : 

The distribution, $2(2fJ^)f in relative units, vs radius, Z, in microns, for 
z =-47m and z = o. 
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Figures H and 12; 

The distribution, $(z,Jl) (solid lines), in relative units, vs radius, £., 
in microns, for various depths, in meters, fron a distributed soxarce, s(z,£) 
(dcished lines; equations 61 and 62) . The same gas diffusion ircdel used to 
calculate $ in Figure 5 is used in this exaitple. 
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BUBBLE RADIUS / (MICRONS) 

Figure 6 
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BUBBLE RADIUS^ (MICRONS) 



■in .. z =- 122 meters 
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BUBBLE RADIUS J (MICRONS) 

Figure 8 
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BUBBLE RADIUS / (MICRONS) 

Figure 9 
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BUBBLE RADIUS Z (MICRONS) 

Figure 10 
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